      PROGRAM LUMM
c     IMPLICIT NONE
      INTEGER N
      REAL A(:,:), lu(:,:), b(:,:)
      REAL COL (:)
cmf$  layout col (:serial)
      INTEGER I,J, K
      real maxdiff, eps
      real rand
      real x

      print *, 'What is the size of the Matrix : '
      read *, N

      allocate (a(1:N,1:N), b(1:N,1:N), lu(1:n,1:n), col(n))

      print *, 'Arrays are allocated'

      call cmf_random (A)

!HPF$ INDEPENDENT, LOCAL_ACCESS
      DO J=1,N 
!HPF$ INDEPENDENT, LOCAL_ACCESS
          DO I=1,N
             lu (i,j) = a(i,j)
          END DO
      END DO 

      PRINT *,'PROGRAM STARTS'
      DO K=1,N-1
         x = lu(k,k)
         do j = k+1, n
            lu(j,K) = lu(j,K)/x
         end do
         col = lu(1:n,k)
!HPF$ INDEPENDENT, LOCAL_ACCESS
         do  j = k+1, n
!HPF$ INDEPENDENT, LOCAL_ACCESS
            do i = k+1, n
              lu (i,j) = lu(i,j) - col(i) * lu(k,j)
            end do
         end do
      END DO

      PRINT *, 'LU ready'
      PRINT *, 'now finds Y with L * Y = E '

!HPF$ INDEPENDENT, LOCAL_ACCESS
      DO J=1,N
!HPF$ INDEPENDENT, LOCAL_ACCESS
          DO I=1,N
              if (I .eq. J) then
                 B(I,J) = 1
               else
                 B(I,J) = 0
              end if 
          END DO
      END DO 

      do j=1, n
         col = lu(1:n,j)
!HPF$ INDEPENDENT, LOCAL_ACCESS
         do k = 1,n
!HPF$ INDEPENDENT, LOCAL_ACCESS
            do i = j+1,n
              b(i,k) = b(i,k) - col(i)*b(j,k)
            end do
         end do
      end do

      PRINT *, 'now finds X with R * X = E '

      do j=n, 1, -1
         col = lu(1:n,j)
!HPF$ INDEPENDENT, LOCAL_ACCESS
         do k = 1,n
            b(j,k) = b(j,k) / col(j)
!HPF$ INDEPENDENT, LOCAL_ACCESS
            do i = 1, j-1
              b(i,k) = b(i,k) - col(i) *b(j,k)
            end do
         end do
      end do

      PRINT *, 'b is the inverse '

c     print *, 'a = ', a
c     print *, 'b = ', b

      PRINT *, 'lu = a * b '

c     lu = 0
!HPF$ INDEPENDENT, LOCAL_ACCESS
      DO J=1,N
!HPF$ INDEPENDENT, LOCAL_ACCESS
          DO I=1,N
             lu (i,j) = 0
          END DO
      END DO 

      do k=1,n
         col = a(1:n,k)
!HPF$ INDEPENDENT, LOCAL_ACCESS
         do j=1,n
!HPF$ INDEPENDENT, LOCAL_ACCESS
           do i=1,n
             lu(i,j) = lu(i,j) + col(i) * b(k,j)
           end do
         end do
      end do
          
      PRINT *, 'matmul correct, now check = E'

c     print *, 'E = ', lu

!HPF$ INDEPENDENT, LOCAL_ACCESS
      do i=1,N
!HPF$ INDEPENDENT, LOCAL_ACCESS
         do j=1,N
            if (i .eq. j) then
               a(i,j) = lu(i,j) - 1
             else
               a(i,j) = lu(i,j)
            end if
            a(i,j) = abs (a(i,j))
         end do
      end do

      maxdiff = 0.0
!HPF$ INDEPENDENT, LOCAL_ACCESS
      do i=1,N
!HPF$ INDEPENDENT, LOCAL_ACCESS
         do j=1,N
           REDUCE (maxval, maxdiff, a(i,j))
         end do
      end do

      deallocate (col, lu, b, a)
      print * , 'maximal eps = ', maxdiff
      END
